/*This script generates the simulation results for the setback illustrations;
Figure 4
Figure A3
Figure A4  

*/


cd "$data"

set more off


*set scalars
global quart_a 160 // sets the acres of the quarter
global quart_l 2640 // sets the length of the square of the quarter
global sec_a 640 // sets the acres of the section
global sec_l 5280 // sets the length of the square of the the section
global padlist 10 // set at 1 acre, but can run over the various pad sizes in acres. 0.06 to .1 is from NREL's Denholm, Hand, Jackson, and Ong (2009)

*make graph
clear
set obs 2000 // sets the number of feet considered for setbacks

gen seqnum=_n

gen set=seqnum/1 

foreach num of numlist $padlist {
	gen pad`num'=`num'/100
	gen d_pad`num'=((2*(((pad`num'*43560)^0.5))^2)^0.5)/2 // converts acres to square feet to calculate the half diagonal length of the pad in feet
	}

gen t_corn_dist= set*2^.5 // sets the distance of the turbine from the corner of the parcel based on the setback (from the sides)

foreach s in quart sec { 
	gen t_cent_dist_`s'= (2*(($`s'_l/2)^2))^(0.5)-t_corn_dist //calculates the distance of the turbine from the center
	
	foreach d of numlist $padlist {
		gen p_cent_dist_`s'_`d'=(t_cent_dist_`s'^2+d_pad`d'^2)^0.5 //sets the distance of the "wide" corners of the pad to the center
		
		gen angle_`s'_`d'=2*asin(d_pad`d'/p_cent_dist_`s'_`d') //sets angle off limits based on the entire diagonal.
		
		local r (($`s'_l)/2) // sets radius of the center pivot

		global a 2 // quadratic formula parameter 
		global b (2*(t_cent_dist_`s'-d_pad`d')) // quadratic formula parameter 
		global c ((t_cent_dist_`s'-d_pad`d')^2-(`r')^2) // quadratic formula parameter 
		

		gen d_`s'_`d'=  (-($b )+(($b )^2-4*(($a )*($c )))^(0.5))/(2*($a )) // finds the distance of the pad cuts into the circle before the full diagonal
		
		replace angle_`s'_`d'=2*asin(d_`s'_`d'/`r') if p_cent_dist_`s'_`d'>`r' //resets the angle when only part of the pad is in the radius of the center pivot

	
		gen acres_`s'_`d'=($`s'_l/2)^2*angle_`s'_`d'/2/43560 //converts the angle to area of circle reduced
		

		gen cost_`s'_`d'=acres_`s'_`d'*2774.76 //price of CPIS land over dryland, estimated by cooley et al. and USDA 2021 value reports

		gen nocropcost_`s'_`d'=acres_`s'_`d'*4624.76 //price of CPIS land estimated by cooley et al. and USDA 2021 value reports
		
		gen perc_`s'_`d'=4*acres_`s'_`d'/($`s'_a*_pi)*100 // finds the percentage of the center pivot taken out of production
		
		gen out_`s'_`d'=(t_cent_dist_`s'-d_pad`d')>$`s'_l/2 // finds whether the pad is entirely beyond the radius of the center pivot
		gen in_`s'_`d'=out_`s'_`d'==0 // inverse of being beyond the center pivot
		
		foreach v in acres cost nocropcost perc {
			replace `v'_`s'_`d'=`v'_`s'_`d'*in_`s'_`d' // makes effect zero for when the pad is beyond the center pivot radius
			replace `v'_`s'_`d'=. if t_cent_dist_`s'-d_pad`d'<0 // makes the effect vanish when the corner of the pad reaches the center

			}
		
		gen s_perc_`s'_`d'=pad`d'/$`s'_a*100 //creates the effect (percent) of the pad on simply irrigated land
		gen s_acres_`s'_`d'=pad`d' //creates the effect (acres) of the pad on simply irrigated land
		gen s_cost_`s'_`d'=pad`d'*1694.80 //price of irrigated land in the states from USDA 2021 report and percentage estimated by cooley et al.

		foreach v in acres cost perc {
				replace s_`v'_`s'_`d'=. if t_cent_dist_`s'<0  //vanishes the effect (percent) of the pad on simply irrigated land once it can't be done
				}	

		}
	}

	

local ytitle_acres "Acres Reduced"
local ytitle_cost "Reservation Price"
local ytitle_perc "Percent of Field"
local ylabel_acres 0(5)30
local ylabel_cost 0(20000)100000
local ylabel_perc 0(5)25
local nam_acres "fig4"
local nam_cost "figA4"
local nam_perc "figA3"
local subnam_quart "a"
local subnam_sec "b"


foreach m in acres cost perc {
	foreach size in  quart sec {

		twoway (line `m'_`size'_10 set, lcolor(forest_green)) (line s_`m'_`size'_10 set, lcolor(blue) lpattern(dash)) if set<1400 ///
		, scheme(s1color) plotregion(lstyle(none)) legend(label(1 "CPIS") label(2 "Irrigation") region(lcolor(white)) ) ytitle(`ytitle_`m'') xtitle("Setback (feet)")  ylabel(`ylabel_`m'') //yscale(range(0(5)50))
			graph export "$results/figures/`nam_`m''`subnam_`size''_`m'_`size'_setback.pdf", replace
			}
		}

/*Extract values for in text references*/
log using "$results/intext/setback_numbers", text replace

	sum acres_quart_10 perc_quart_10 cost_quart_10 nocropcost_quart_10 if set == 390
	sum acres_quart_10 perc_quart_10 cost_quart_10 nocropcost_quart_10 if set == 1000
	sum acres_quart_10 perc_quart_10 cost_quart_10 nocropcost_quart_10 if set == 1250
	sum acres_quart_10 perc_quart_10 cost_quart_10 nocropcost_quart_10 if set == 1285

log close


